FAST ALGORITHMS FOR COMPUTING THE BOLTZMANN 

COLLISION OPERATOR 



CLEMENT MOUHOT AND LORENZO PARESCHI 



Abstract. The development of accurate and fast numerical schemes for the five fold 
Boltzmann collision integral represents a challenging problem in scientific computing. 
For a particular class of interactions, including the so-called hard spheres model in di- 
mension three, we are able to derive spectral methods that can be evaluated through fast 
algorithms. These algorithms are based on a suitable representation and approximation 
of the collision operator. Explicit expressions for the errors in the schemes are given and 
spectral accuracy is proved. Parallelization properties and adaptivity of the algorithms 
are also discussed. 
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1. Introduction 

The Boltzmann equation describes the behavior of a dilute gas of particles when the 
only interactions taken into account are binary elastic collisions. It reads for x,v £ M. d 
(d>2) 

^ + v-V x f = Q(f,f) 

where f(t,x,v) is the time-dependent particle distribution function in the phase space. 
The Boltzmann collision operator Q is a quadratic operator local in (t,x). The time and 
position acts only as parameters in Q and therefore will be omitted in its description 

(1.1) Q(fJ)(v)= [ B(\v-v.\,cob9) [fif'-Uf) dv*dcr. 
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In we used the shorthand / = f(v), /* = /(t>*), /' = f(v'), f # = f(v'^). The velocities 
of the colliding pairs (v,v*) and (v' ,v' # ) are related by 

, v + v* \v — vj , v + v* \v — vJ 

v = 1 cr, = a. 

2 2 ' * 2 2 

The collision kernel B is a non-negative function which by physical arguments of invariance 
only depends on \v — v*\ and cos# = g ■ a (where g = (v — v*)/\v — v*\). 

Boltzmann's collision operator has the fundamental properties of conserving mass, mo- 
mentum and energy 



Q(f, f)4>{v) dv = 0, (j)(v) = l,v 1 , . . . ,v d ,\v 

and satisfies the well-known Boltzmann's H theorem 

d 



2 



dt 



[ f log fdv = -[ Q(/,/)log(/)d«>0. 



The functional — J f log / is the entropy of the solution. Boltzmann's H theorem implies 
that any equilibrium distribution function, i.e. any function which is a maximum of the 
entropy, has the form of a locally Maxwellian distribution 

p ( \u — v 1^ 

M(p '"' T)(v) = (^rj^ exp \ — 

where p, u, T are the density, mean velocity and temperature of the gas, defined by 
p = [ f(v)dv, u = - [ vf(v)dv, T = — f \u — v\ 2 f(v) dv. 

jR d P JR d dp J R d 

For further details on the physical background and derivation of the Boltzmann equation 
we refer to [T5| linj. 

The construction of numerical methods for Boltzmann equations represents a real chal- 
lenge for scientific computing and it is of paramount importance in many applications, 
ranging from rarefied gas dynamics (RGD) plasma physics [IH1, granular flows [2112]) 
semiconductors [2131 an d quantum kinetic theory |27Tj . 

Most of the difficulties are due to the multidimensional structure of the collisional 
integral Q, as the integration runs on a 5-dimensional unflat manifold. In addition to the 
unpracticable computational cost of deterministic quadrature rules the integration has to 
be handled carefully since it is at the basis of the macroscopic properties of the equation. 
Additional difficulties are represented by the stiffness induced by the presence of small 
scales, like the case of small mean free path |26j or the case of large velocities |22| . 

For such reasons realistic numerical computations are based on probabilistic Monte- 
Carlo techniques at different levels. The most famous examples are the direct simulation 
Monte Carlo (DSMC) methods by Bird 0] and by Nanbu |35j . These methods preserve 
the conservation properties of the equation in a natural way and avoid the computational 
complexity of a deterministic approach. However avoiding the low accuracy and the fluc- 
tuations of the results becomes extremely expensive in presence of nonstationary flows or 
close to continuum regimes. 
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Among deterministic approximations, one of the most popular methods in RGD is 
represented by the discrete velocity models (DVM) of the Boltzmann equation. These 
methods ^H^0E1EH1I1Z1 are based on a regular grid in the velocity field and construct 
a discrete collision mechanics on the points of the grid in order to preserve the main 
physical properties. Unfortunately DVM have the same computational cost of a product 
quadrature rule and due to the particular choice of the nodes imposed by the conservation 
properties the accuracy of the schemes seems to be less than first order |37l 1361 05] . 

More recently a new class of methods based on the use of spectral techniques in the ve- 
locity space has attracted the attention of the scientific community. The method was first 
developed for kinetic equations in [10], inspired from spectral methods in fluid mechan- 
ics and the use of Fourier transform tools in the analysis of the Boltzmann equation 
It is based on a Fourier-Galerkin approximation of the equation. Generalizations of 
the method and spectral accuracy have been given in |411 142j . This method, thanks to 
its generality, has been applied also to non homogeneous situations [21], to the Landau 
equation |221 I43| and to the case of granular gases |34l I23j . A related numerical strategy 
based on the direct use of the fast Fourier transform (FFT) has been developed in JS][H]- 

The lack of discrete conservations in the spectral scheme (mass is preserved, whereas 
momentum and energy are approximated with spectral accuracy) is compensated by its 
higher accuracy and efficiency. In fact it has been shown that these spectral schemes permit 
to obtain spectrally accurate solutions with a reduction of the computational cost strictly 
related to the particular structure of the collision operator. A reduction from 0(N 2 ) to 
0(N log 2 N) is readily deducible for the Landau equation, whereas in the Boltzmann case 
such a reduction had been obtained until now only at the price of a poor accuracy (in 
particular the loss of the spectral accuracy), see [HUB]- 

Finally we mention that spectral methods have been successfully applied also to the 
study of non cut-off Boltzmann equations, like for RGD in the grazing collision limit |45| 
and for granular flows in the quasi-elastic limit (31]. In particular, during these asymptotic 
processes it is possible to obtain intermediate approximations that can be evaluated with 
fast algorithms that brings the overall computational cost to 0(N log 2 N). These idea has 
been used in |39j to obtain fast approximated algorithms for the Boltzmann equation. 

For a recent introduction to numerical methods for the Boltzmann equation and related 
kinetic equations we refer the reader to |19j . 

In this paper we shall focus on the two main questions in the approximation the Boltz- 
mann equation by deterministic schemes, that is the computational complexity and the 
accuracy of the numerical schemes for computing the collision operator Q. 

Let us mention that a major problem associated with deterministic methods that use 
a fixed discretization in the velocity domain is that the velocity space is approximated 
by a finite region. Physically the domain for the velocity is M. d . But, as soon as d > 2, 
the property of having compact support is not conserved by the collision operator (in 
fact for some Boltzmann models in dimension d = 1, like granular models, the support 
is conserved |34| ) . In general the collision process "spreads" the support by a factor v2 
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(see |461 1321 ). As a consequence, for the continuous equation in time, the function / is 
immediately positive in the whole velocity domain Mr. 

Thus at the numerical level some non physical condition has to be imposed to keep the 
support of the function in velocity uniformly bounded. In order to do this there are two 
main strategies, which we shall make more precise in the sequel. 

(1) One can remove the physical binary collisions that will lead outside the bounded 
velocity domain, which means a possible increase of the number of local invari- 
ants. If this is done properly (i.e. "without removing too many collisions"), the 
scheme remains conservative (and without spurious invariants). However this trun- 
cation breaks down the convolution-like structure of the collision operator, which 
requires the invariance in velocity. Indeed the modified collision kernel depends on 
v through the boundary conditions. This truncation is the starting point of most 
schemes based on discrete velocity models in a bounded domain. 

(2) One can add some non physical binary collisions by periodizing the function and the 
collision operator. This implies the loss of some local invariants (some non physical 
collisions are added). Thus the scheme is not conservative anymore, except for the 
mass if the periodization is done carefully (and possibly the momentum if some 
symmetry properties are satisfied by the function). In this way the structural 
properties of the collision operator are maintained and thus they can be exploited 
to derive fast algorithms. This periodization is the basis of the spectral method. 

Note that in both cases by enlarging enough the computational domain the number of 
removed or added collisions can be made negligible (as it is usually done for removing the 
aliasing error of the FFT, for instance see JJ) as well as the error in the local invariants. 

In this paper we shall focus on the second approach, which means that the schemes 
have to deal with some aliasing error introduced by the periodization. In this way, for a 
particular class of interactions, using a Carleman-like representation of the collision oper- 
ator we are able to derive spectral methods that can be evaluated through fast algorithms. 
The class of interactions includes Maxwellian molecules in dimension two and hard spheres 
molecules in dimension three. 

The rest of the paper is organized in the following way. In Section [2] we introduce a 
Carleman-like representation of the collision operator which is used as a starting point for 
the development of our methods. After the derivation of the schemes the details of the 
fast spectral algorithm together with its accuracy properties are given in Sectional In a 
separate Appendix we show a possible way to extend the present fast schemes to general 
collision interactions. 

2. Carleman-like representation and approximation of the collision 

operator 

In this section we shall approximate the collision operator starting from a representation 
which somehow conserves more symmetries of the collision operator when one truncates 
it in a bounded domain. This representation was used in [SI |SJ |51 and it is close to the 
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classical Carleman representation (cf. ^21)- Also the kind of periodization inspired from 
this representation was implicitly used in |S]. 

2.1. The Boltzmann collision operator in bounded domains. The basic identity 
we shall need is 

(2.1) - / F(\u\a - u) do- = , , , „ / 5(2x ■ u + Id 2 ) F(x) dx, 

2 J§d-i \u\ d 1 J Rd 

and can be verified easily by completing the square in the delta Dirac function, taking the 
spherical coordinate x = r a and performing the change of variable r 2 = s. 
Setting u = v — v * we can write the collision operator in the form 



Q(f,f)(v) 



B(|u|,cos0) 



/K-^)/f» + ^)-/M/M 



and thus equation (|2.1j) yields 
Q(fJ)(v) = 2 [ \f B (\u\,^\ -^5(2x.u + \xf) 



da > dv* 



f(v.-x/2)f(y + x/2)-f(v.)f(y) 
Now let us make the change of variable x — > x/2 in x to get 



dx > dv*. 



Q(f,f)(v) = 2 



d+l 



v,eM d Jx 



B \u 



x ■ u 



1 



Id— 2 



5(4 x ■ u + 4|a?| 



|x||n|/ \u\ 

[/(«* - x) f(v + x) - f(v*) f(v)] dxdv* 

and then setting y = — v — x in v* we obtain 



Q(fJ)(v) = 2 



d+l 



B \u 



x ■ u 



1 



Id— 2 



5(—4x ■ y) 



\X\\U\ J \u\ 

[f(v + y) f(v + x) - f(v + x + y) f(v)] dx dy 
where now u = — (x + y) . Thus in the end we have 



Q(f,f)(v) = 2 



d-l 



1 


/ H 







1 



x ■ (x + y) 



\x\ \x + y\ J \x + y 



Id— 2 



5(x -y) [f(v + y) f(v + x) - f(v + x + y) f(v)] dx dy. 



Figure n sums up the different geometrical quantities of the usual representation and the 
one we derived from Carleman 's one. 
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Figure 1. Geometry of the collision <-> (i/,t>*). 



Now let us consider the bounded domain T>t = [—T,T] d (0 < T < +oo). There are 
two possibilities of truncation to reduce the collision process in a box. From now on let 
us write 

B(x, y) = I*' 1 b(\x + V\- X ',\ X + V 1 ) \x + y\- {d - 2) . 
V \x\\x + y\J 

One can easily see that on the manifold defined by x ■ y = 0, a simpler formula is (using 
the parities of the collision kernel) 

(2.2) B{x,y) = B{\x\,\y\) = 2 d ~ l B ( y/\x\* + \y\\ }f ) (|x| 2 + |y| 2 )-^. 

V v\ x \ +\y\ J 

First one can remove the collisions connecting with some points out of the box. This is the 
natural preliminary stage for deriving conservative schemes based on the discretization of 
the velocity. In this case there is no need for a truncation on the modulus of x and y since 
we impose them to stay in the box. It yields 

Q tT (f,f)(v)= I I 1 B(x,y)5(x-y) 

J J{x,y£R d | v+x,v+y,v+x+y <=T> T j 

[f(v + y) f(v + x) - f(v + x + y) f(v)] dx dy 
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defined for v £ T>t- One can easily check that the following weak form is satisfied by this 
operator 

(2.3) ( Q tl {fJ)^{v)dv = - III B(x,y)5(x-y) 

J 4 J J J {v,x,y£R d | v,v+x,v+y,v+x+y GV T j 

f(v + x + y) f(v) [tp(v + y) + (p(v + x) — ip(v + x + y) — <p(v)] dv dx dy 

and this implies conservation of mass, momentum and energy as well as the H theorem 
on the entropy. Note that at this level this formulation gives no advantage with respect 
to the usual one obtained from (jl.ljl by restricting v, «*, v', v'* £ T>t (except that consis- 
tency results for discrete velocity models seem easier to prove when they are derived by 
quadrature on this formulation, see [SHI)- The problem of this truncation on a bounded 
domain is the fact that we have changed the collision kernel itself by adding some artificial 
dependence on v, v*, v', u*. In this way convolution-like properties are broken. 

A different approach consists in periodizing the function / on the domain T>x- This 
amounts in adding some non-physical collisions by connecting some points in the domain 
T>t which are geometrically included in a collision circle "modulo T" (i.e. up to a transla- 
tion of T of certain points in certain directions). Here we have to truncate the integration 
in x and y since periodization would yield infinite result if not. Thus we set them to vary 
in Br, the ball of center and radius R. For a compactly supported function / with 
support Bs, we take R = 2S in order to obtain all possible collisions. Then a geometrical 
argument (see jl^) shows that using the periodicity of the function it is enough to take 
T > (1 + 3\/2)5*/2 to prevent intersections of the regions where / is different from zero. 
Note that here this so-called dealiasing condition is slightly worst from the one in |41j . 
since the truncation on the modulus of x and y in the ball Br implies only a truncation 
in the ball B^ R for the relative velocity. 

The operator now reads 

(2.4) Q R (f,f)(v)= [ [ B(x,y)S(x-y) 

JxgBr Jy&B R 

[f(v + y)f(v + x) - f(v + x + y)f(v)] dx dy 

for v G T>t (the expression for v G M rf is deduced by periodization). The interest of this 
representation is to preserve the real collision kernel and its properties. 

By making some translation changes of variable on v (by x, y and x + y), using the 
changes x — > — x and y — > —y and the fact that 

B(-x, y) 5{-x ■ y) = B(x, y) 5(x ■ y) = B(x, -y) 5(x ■ -y) 

one can easily prove that for any function <p periodic on T>t the following weak form is 
satisfied 

(2.5) / Q R {f,f)ip(v)dv = \ III B(x,y)5{x-y) 
JV T 4 Jv&Vt Jx£B r Jy<=B R 

f(v + x + y)f(v) [ip(v + y) + ip(v 4- x) — <p(v + x + y) - (f(v)] dv dx dy. 
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About the conservation properties one can shows that 

(1) The only invariant tp is 1: it is the only periodic function on T>t such that 

(p(v + y) + ip(v + x) — ip(v + x + y) — ip(v) = 

for any v E T>t and X-Ly E Br (see |13| for instance). It means that the mass is 
locally conserved but not necessarily the momentum and energy. 

(2) When / is even there is global conservation of momentum, which is in this case. 
Indeed Q R preserves the parity property of the solution, which can be checked 
using the change of variable x ^ —x, y ^ —y. 

(3) The collision operator satisfies formally the H theorem 

/ Q R (f,f)\og(f)dv<0. 

(4) If / has compact support included in Bg, and we have R = 2S and T > (3\/2 + 
1)572 (no aliasing condition, see |41j for a detailed discussion), then no unphysical 
collisions occur and thus mass, momentum and energy are preserved. Obviously 
this compactness is not preserved with time since the collision operator spreads 
the support of / by a factor y/2. 

To sum up one could say that the lack of conservations originates from the fact that the 
geometry of the collision does not respect the periodization. 

Finally we give the Cauchy theorems for the homogeneous Boltzmann equations in T>t 
computed with Q^ r or Q R . 

Theorem 2.1. Let /o G L 1 (2?t) be a nonnegative function. Then there exists a unique 
solution f G C 1 (M+, L 1 (I?t)) to the Cauchy problems 

(2.6) % = Qtr V>fi> /(* = M=/o 

(2-7) ^ = Q R (f,f), /(t = 0,.)=/o 

which is nonnegative and has constant mass (and so constant L 1 norm). If fo has finite 
entropy, the entropy is finite and non- decreasing for all time. Moreover in the case (|2.(i|) . if 
/o has finite momentum (respectively energy) onT>T, the momentum (respectively energy) 
is conserved with time. 

Remark: When the initial data fo is nonnegative and has finite mass and entropy, it is 
possible to show by the Dunford- Pettis compactness theorem that the solution / converges 
weakly in L 1 (Pj>), as t goes to infinity, to the unique maximum of the entropy functional 
compatible with the conservation law(s) (and the periodicity in the case (|2.7|) ). In the 
case (|2.6|) this equilibrium state is a sort of truncated Maxwellian on T>t defined by the 
conservation laws (see |H3)- In the case (|2.7|) this equilibrium state is a constant defined 
by the mass of the initial data, which is due to the effect of aliasing in the very long-time. 
We omit the proof for brevity. 
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Proof of Theorem \2.1\ For clarity we briefly sketch the main lines of the proof. The ex- 
istence and uniqueness are proved by the method of Arkeryd for bounded collision ker- 
nels, see |U Part I, Proposition 1.1]. In our case the collision kernel is bounded because 
of the boundedness of the domain. The only a priori estimate required in ^ Part I, 
Proposition 1.1] is the mass conservation, valid for the two equations under consideration. 
This method is based on a monotonicity argument to prove propagation of the sign of 
the solution. The argument relies on a splitting of the collision operator Q into a gain 
part Q + which is monotonic {i.e. Q + (/, /) is non-negative when / is non-negative), and 
a loss part Q~ which writes Q~(f,f) = L(f)f with L is a linear operator such that 
||^(/)||oo < C II/Uli- One can check easily that this splitting is still valid for the two colli- 
sion operators Q^ T and Q R . For brevity we omit the details and refer to the article pp. The 
conservation law(s) and the H theorem are deduced from the weak forms (|2.3|) and (|2.5|) 
(see the proof of ^ Part I, Proposition 1.2] and ^ Part I, Theorem 2.1]). □ 

2.2. Application to spectral methods. In this Section we use the representation Q R 
to derive new spectral methods. The spectral methods for kinetic equations originated 
in the works of |4()1 141 j , and were further developed in |421 12*1] . Before they had a long 
history in fluid mechanics, see jllj . 

The main change compared to the usual spectral method is in the way we truncate the 
collision operator. In fact as we shall see in the next section this yields better decoupling 
properties between the arguments of the operator. 

To simplify notations let us take T = ir. Hereafter we use just one index to denote the 
d-dimensional sums of integers. 

The approximate function f^ is represented as the truncated Fourier series 

N 

fN(v) = f* eik ' V > 



k=-N 



1 

{2n) c 



-ik-v 

J 

!p.. 

The spectral equation is the projection of the collision equation in F N , the (2N + l) d - 
dimensional vector space of trigonometric polynomials of degree at most N in each direc- 
tion, i.e. 

df N 



V N Q H (f N ,f. 



N 



dt 

where Vn denotes the orthogonal projection on in L 2 ('D 7r ). A straightforward compu- 
tation leads to the following set of ordinary differential equations on the Fourier coefficients 



(2.8) f' k (t)= Yl PMflL, k = -N,...,N 

l,m=-N 

l-\-m—k 
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where (3(1, m) are the so-called kernel modes, given by 



0(1, m) 



<x&B R Jy&B R 
The kernel modes can be written as 



B(x,y)5(x ■ y) 



£ il-x e im-y _ e im-(x+y) 



dx dy. 



(3(1, m) = (3(1, m) — (3(m, m) 

where 

(3(1, m)= ! [ B(x,y)5(x-y)e il - x e im - y dxdy. 
JxeB R JyeB R 

Therefore in the sequel we shall focus on (3, and one easily checks that (3(1, m) depends 
only on \m\ and |Z • m\. 

Note that the usual way to truncate the Boltzmann collision operator for periodic 
function starts from the following representation (see |41| ) 



da du 



(2.9) Q(f,f)= I I B(\u\,co S e) 



f(v -(u- \u\a)/2)f(v -(u + \u\a)/2) - f(v)f(v - u) 
and then truncate the parameter u = x + y in order that u E Br. Thus we have 

QmnMJ)( v )= / B(x,y)5(x ■y)x{\ x+y \<R} 

Jx£R d Jy£R d 

[f(v + y)f(v + x) - f(v + x + y)f(v)] dx dy 

where X{\x+y\<R} denotes the characteristic function of the set {\x + y\ < R}. One can 
notice that here x and y are also restricted to the ball Br but the condition \x + y\ 2 = 
\ x \ 2 + 1 2/ 1 2 ^ R 2 couples the two modulus, such that the ball is not completely covered (for 
instance, if x and y have both modulus R, the condition is not satisfied, since \x + y\ = 
y/2R). 

Finally let us compare the new kernel modes with the usual ones. As a consequence of 
the representation (|2.9j) . the usual kernel modes (cf. are 



B(\u\,cos6) 



■ u-(Z+m) + |u|cr-(m — I) 



-i(u-m) 



'u£B R Ja& d - 1 

and hence coming back to the representation in x and y, 



/3usual(^"l) 



x£B R Jy£B R 



B(x,y)8(x ■ y)X{\x+y\<R} 



e il-x e %m-y _ e im-(x+y) 



da du 



dx dy. 



Thus the usual representation contains more coupling between x and y and it is less 
appropriate for the construction of fast algorithms. 
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3. Fast spectral algorithm for a class of collision kernels 

As soon as one is searching for fast deterministic algorithms for the collision operator, 
i.e. algorithm with a cost lower than 0(N 2d+£ ) (which is the cost of a usual discrete 
velocity model, with typically e = 1), one has to find some way to compute the collision 
operator without going through all the couples of collision points during the computation. 
This leads naturally to search for some convolution structure (discrete or continuous) in the 
operator. Unfortunately, as discussed in the previous sections, this is rather contradictory 
with the search for a conservative scheme in a bounded domain, since the boundary con- 
dition needed to prevent for the outgoing or ingoing collisions breaks the invariance. Thus 
fast algorithms seem more adapted to spectral methods, or more in general to methods 
where the invariance is conserved thanks to the periodization. 

Here we search for a convolution structure in the equations (|2.8|) . The aim is to ap- 
proximate each (3(1, m) by a sum 

A 

0(1, m) ^^2a p (l)a p (m). 
P =i 

This gives a sum of A discrete convolutions and so the algorithm can be computed in 
0(A N d log 2 N) operations by means of standard FFT techniques |11[ I16|. Obviously this 
is equivalent to obtain such a decomposition on (3. To this purpose we shall use a further 
approximated collision operator where the number of possible directions of collision is 
reduced to a finite set. 

The starting point of our study is an idea of |Hj: use the Carleman-like representation 
(|2.4|) to obtain a convolution structure for every fixed directions of the vectors x and y. 
In this work jg] the corresponding set of directions 

S = {(e,e')€§ N ~ 1 xS N - 1 | e±e'} 

is very difficult to discretize in a way that preserves the symmetry properties of the collision 
operator. No systematic process is available and the discretization is done only for some 
particular number of grid points. Then the FFT is used in each couple of direction and 
finally a correction is imposed at the end to preserve the conservation laws. However no 
consistency result is available and the accuracy suggested by the numerical simulations is 
of order 1. The two main new ingredients of our method are: 

• First we project the collision operator on the Fourier basis. This enables to inte- 
grate one of the two coordinates of the manifold S and to reduce to the discretiza- 
tion of the sphere S^ -1 . This discretization is straightforward and can be made 
easily to preserve the symmetries of the collision operator. Moreover it reduces the 
complexity of the algorithm by suppressing N — 2 degrees of freedom to discretize. 

• Second we choose to discretize S^ -1 by the rectangular rule. Indeed the periodiza- 
tion shall imply that this quadrature rule is of infinite order. This point will allow 
to obtain a spectrally accurate scheme, and adaptativity properties. 
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3.1. A semi-discrete collision operator. We write x and y in spherical coordinates 
(3.1) Q R (f,f)(v) = ] [ [ S(e-e')dede' 

1 fa f R ^ {p,)d 2 ^ + ^^^^ + ^ " + P& + dRdP \ 



Let us take A a set of orthogonal couples of unit vectors (e, e'), which is even: (e, e') £ A 
implies that (— e, e'), (e, — e') and (— e, — e') belong to „4 (this property on the set „4 is 
required to preserve the conservation properties of the operator). Now we define to 
be 



1 



4 J(e,e')eA 



R r R 



-R J-R 



[f(v + p'e')f(v + pe) - f(v + pe + p'e')f(v)] dp dp' \ dA 



where dA denotes a measure on A which is also even in the sense that dA(e, e') = 
dA(—e,e') = dA(e,—e') = dA(—e,—e'). Using translation changes of variable on v by 
pe, p'e' and pe + p'e' and the symmetries of the set A one can easily derive the following 
weak form on Q^. For any function (p periodic on T>t, 

l Q R ' A (fJ)<p(v)dv=±; [ [ [ R [ R p d - 2 (p') d - 2 B(p, P ') 

JV T iD Jv£V T J(e,e')eAJ-RJ-R 

f(v + pe + p'e')f(v) ip(v + p'e') + (p(v + pe) — ip(v + pe + p'e') — ip(v) dpdp' dAdv. 

This immediately gives the same conservations properties as Q R . Of course one could also 
prove exactly as for Q R : 

Theorem 3.1. Let G L 1 (T>t) be a nonnegative function. Then there exists a unique 
solution f G C 1 (R+, L 1 (Dt)) to the Cauchy problem 

% = Q R ' A (f,f), f(t = o,-) = fo 

which is nonnegative and has constant mass (and so constant L 1 norm). Moreover, if fo 
has finite entropy, the entropy is non- decreasing with time. 

3.2. Expansion of the kernel modes. We make the decoupling assumption that 

(3.2) B(x,y)=a(\x\)b(\y\). 

This assumption is obviously satisfied if B is constant. This is the case of Maxwellian 
molecules in dimension two, and hard spheres in dimension three (the most relevant kernel 
for applications). Extensions to more general interactions are discussed in the Appendix. 
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First let us deal with dimension 2 with B = 1 to explain the method. Here we write x 
and y in spherical coordinates x = pe and y = p'e' to get 



(3(1, m) = 
Let us denote by 



eGS 1 Je'eS 1 



5(e ■ e') 



R 



R 



e ip{ - le) dp 



R 



R 



Jp>(m-e>) df)l 



de de . 



R 



e ips dp, 



for s£K. It is easy to see that (p R is even and we can give the explicit formula 



4> 2 R (s) = 2RSmc(Rs) 



with Sinc(0) = (sin0)/0. 
Thus we have 

0(1, m) 



eeS 1 Je'eS 1 
2 



S(e ■ e') (f> R (l ■ e) 4> 2 R (m ■ e') dede' 



and thanks to the parity property of (j) R we can adopt the following periodic parametriza- 
tion 

0(1, m) = <f) R (l ■ e e ) 4> R (m • e e+7r / 2 ) dO. 
Jo 

The function 8 — > 4> R (l ■ e$) 4>R( m " e e+ir/2) is periodic on [0,tt] and thus the rectangular 
quadrature rule is of infinite order and optimal. A regular discretization of M equally 
spaced points thus gives 

M-l 

Khm) = — a p (l)a' p (m) 

p=0 

with 

«p(0 = 4>r(1 ■ %)> «p( m ) = 4> R (m ■ e 9p+7r/2 ) 

where 8 p = irp/M. 

More generally under the decoupling assumption (jS2J) on B, we get the following de- 
composition formula 

M-l 

0(1, m) = — a p (l)a p (m) 

p=0 



where 
and 



MO = (t>R,a{ 1 ■ %)> a ' P ( m ) = <l>R,b( m ■ e e p +n/2) 



<t>R,a( s ) 

with Op = irp/M. 



R 



a(p)e ips dp, 4>l b (s) 



R 



b(p') e ip s dp' 
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Remark: In the symmetric case a = b (for instance for hard spheres) it is possible to 
parametrize /3(l,m) as 

/•tt/2 

0(l,m)=2 / 4>% a (l ■ e e ) <f> 2 R J™ ■ ee+ir/2) dO 
Jo 

and the function 9 — > (j>\ a (l-eo) 4>Ra( m ' e d+ir/2) is periodic on [0,7r/2]. Thus the decompo- 
sition can be obtained by applying the rectangular rule on this interval. At the numerical 
level it yields a reduction of the cost by a factor 2. 

Now let us deal with dimension d = 3 with B satisfying the decoupling assumption ()3.2|) . 
First we change to the spherical coordinates 



0(l,m) 



e'en 2 



5{e-e') 



R 



R 



\p\a(p)e i ^dp 



R 



-R 



IpXpV^™- 60 dp' 



de de 



and then we integrate first e' on the intersection of the unit sphere with the plane e -1 -, 

1 



p(l,m) 



ees 2 



e'eS 2 ne-L 



4> R b (m ■ e') de' 



de 



where 

= I \P\ «(P) ^ S dp, 4>r,M = [ |p| b(p) e^ s dp. 

J-R J-R 

Thus we get the following decoupling formula with two degrees of freedom 



0(l,m) 

where denotes the half-sphere and 

V>!, 6 (n e ±(m)) = 



<l>R,a( l - e ) ^fl,b(n e ±(m))de 



HUM! cose) de, 



(this formula can be derived performing the change of variable de' = sin (9 de dip with the 
basis (e, u = II e x (m)/|II e x (m)\, e x u)). 

Again in the particular case where B = 1 (hard spheres model), we can compute 
explicitly the functions 4> R (in this case a = b = 1), 

0|( s ) = R 2 [2Sinc(fls) - Sinc 2 (i?s/2)] . 

Now the function e — ► (j) R a (l ■ e) ip R b (ll e x (m)) is periodic on and so the rectangular 
rule is of infinite order and optimal. Taking a spherical parametrization (0,(p) of e G 
and uniform grids of respective size M\ and M2 for and we get 

2 Mi,M 2 



MiM 2 



where 



«p,?(0 = <f>R,a(l ■ H^)) 



-{9p,<pq) 



m 
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and 

From now on we shall consider this expansion with M = M\ = M2 to avoid anisotropy in 
the computational grid. 

Remarks: 

1. It is possible to give more general exact formula in dimension 2 and 3 when a(r) = |r|*, 
b(r) = |r|* with t, t' G N by computing derivatives along along s of the two quantities 

r-R r-R 
I sm(ps)dp, / cos(ps)dp. 
Jo Jo 



2. For any dimension, we can construct as above an approximated collision operator 
qR,Am w j t h 

A M = {(e,e') G S**" 1 x S**" 1 | e G S^J., e' G e^S^ 1 } 

where S^^_ denotes a uniform angular discretization of the half sphere with M points in 
each angular coordinate (the other half sphere is obtained by parity). Let us remark that 
this discretization contains exactly M d ~ x points. From now on we shall denote 

qR,M = qR,A m = ^ Qp' M . 

p=l 



3.3. Spectral accuracy. In this paragraph we are interested in computing the accuracy 
of the scheme according to the three parameters N (the number of modes), R (the trun- 
cation parameter), and M (the number of angular directions for each angular coordinate). 
Instead of looking at the error on each kernel mode it is more convenient to look at the 
error on the global operator. Here the Lebesgue spaces L p , p = 1 ... + 00, and the periodic 
Sobolev spaces Hp, k = . . . + 00 refer to XV 

In order to give a consistency result, the first step will be to prove a consistency result 
for the approximation of Q R by Q R,M . 

Lemma 3.2. The error on the approximation of the collision operator is spectrally small, 
i.e. for all k > d — 1 such that f G Hp 

\m s j)-Q^( S j)\\,<c f MH ^ H ' . 



Hi 
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Proof of Lemma \3. 6 A Starting from (j3.1|) . one gets 



Q R (gJ)(v) 



R r R 



eeSi" 1 Je'&S d - 1 ne ± J-R J-R 

[g(v + p'e')f{v + pe) - g(y + pe + p'e')f{v)] dp dp' de! 



de. 



As the function in the brackets is a periodic function of e on S^. -1 with period ir in each 
coordinate, one can apply the error estimate for the rectangular rule (see for instance |48| 
Theorem 19.10]). This error estimate is valid for k > d — 1 and depends on the derivative 
along e of this functional on the following way 



\\Q R (gJ)-Q H > M (gJ)\\L* < 



R,M i 



c 



2 k M k 



d-l 

E 

8=1 







e / 6 §d-i ne ± j_ R j_ R 



R r R 



2 / J\d-2 



P d - Z on 



B(p,p') [g(v + p'e')f(v + pe) - g(v + pe + p'e')f(v)] dp dp' de' 



de 



LI 



where the constant is independent on k and d k , is the derivative of order k along the 
coordinate ej. Then a straightforward computation gives 



\\Q R (9,f)-Q H ' M (gJ)h2 < 



R,M / 



CR 



2 k M k 



k d -! 

E 



i=i 



£ (k')(\\Q R ' + (\d k '9\Ad k "f\)\\ L2 

k'+k"=k 

+ \\Q R '-(\d k 'g\,\d k "f\)\\ 



L 2 , 



where d k ' and d k " denote some derivatives of order k! and k" . Then using the estimates 

\\Q R ' + (gJ), Q R '-(9,f)h*<C\\g\\ L 4f\\ L 2 

proved in US! 1 , we get 

C F? k 

\\Q R (g, f) - Q R > M (g, f)\\ L > < \\g\\u- 11/11** 

which concludes the proof. □ 



For the second step we shall use the consistency result |41i Corollary 5.4] on the operator 
Q R , which we quote here for the sake of clarity. 

Lemma 3.3. For all k £ N such that f G , 

\\Q R (f,f)-V N Q R (f N J N )\\ L2 < ^ (||/||^ + \\Q R (fN,f N )\\ H k). 



1 Which are consequences of the L p estimates proved in |27l I28| . and revisited in |33|. 
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Combining these two results, one gets the following consistency result 
Theorem 3.4. For all k > d - 1 such that f E H^V^), 



J- L> \\J1\\ \ Tjk ( ' / 

\\Q R (f,f)-r N Q R ' M (f N ,f N )\\ L , <Cl m P +j±(\\f\\H* + \\Q R (fN,fN)\\H* 




Proof of Theorem \3.4\ By triangular inequality 




\\V N {Q R (f N ,f N )-Q R > M (f N ,f N ))\\ L 2 < \\Q R (f N jN)-Q R ' M (fN,fN)\\L* 




The second term in the right-hand side is controlled by Lemma 13.31 which concludes the 



Now let us focus briefly on the macroscopic quantities. In fact here no additional error 
(related to M) occurs, compared with the usual spectral method, since the approximation 
of the collision operator that we are using is still conservative. First with Lemma 13.21 at 
hand one can establish the estimate 



for a constant uniform in M. Then following the method of [411 Remark 5.4] and using 
this estimate we obtain the following spectral accuracy result 



| (Q R ' M (f, /), <P)-(VnQ R,M Un, fN),<p)\& < prlMks (\\f\\ H ^ + \\Q R ' M (fNjN)\\H*) 



where <p can be replaced by v, \v\ 2 . Indeed there is no need to compare the momenta 
of VnQ R ' M '(/jv, In) with those of Q R (f,f) since Q R,M is also conservative, and so they 
can be compared directly to those of Q R,M . Thus the error on momentum and energy is 
independent on M and is spectrally small according to N even for very small value of the 
parameter M. 

3.4. Implementation of the algorithm. The final spectral scheme depends on the 
three parameters N, R, and M. The only conditions on these parameters is the no- 
aliasing condition that relates R and the size of the box T (here it). A detailed study of 
the influence of the choices of N and R has been done in jlj. Here we are interested only 
in the influence of M over the computations, since M controls the computations speed-up. 



proof. 



□ 



\\Q R > M (gJ)\\ L i<c\\g\\ H a\\f\\ H «, 
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N 


M=2 


M=4 


M=8 


M=16 


32 


2.129E-4 


1.993E-05 


2.153E-05 


2.262E-5 


64 


2.109E-4 


7.122E-10 


6.830E-10 


6.843E-10 


128 


2.112E-4 


3.116E-12 


3.117E-12 


3.117E-12 



Table 1. Relative L\ norm of the error for different values of N and M 
for the fast spectral method. 



The method of the previous subsections yields a decomposition of the collision operator, 
which after projection on ¥ N gives the following decomposition 

(3.3) V N Q R ' M = £ V N Qf M . 

P =i 

Each VnQp' M can be computed with a cost 0(N d log 2 N). Thus for a general choice of 
M and N we obtain the cost 0(M d ~ 1 N d log 2 N). The decomposition ()3.3|) is completely 
parallelizable and thus the cost can be strongly reduced on a parallel machine (theoritically 
up to 0(N d log 2 N)). One just has to make independent computations for the M d ~ l terms 
of the decomposition. 

Moreover the formula of decomposition is naturally adaptive (that is the number M 
can be made space dependent), which can be quite useful in the inhomogeneous setting, 
where some regions deserve less accuracy than others. Since it relies on the rectangular 
formula, whose adaptivity property is well known, one can easily double the number of 
directions M if needed, without computing again those points already computed. 

Finally the decomposition can be also interesting from the storage viewpoint, as the 
classical spectral method requires the storage of a N d x N d matrix whereas our method 
requires the storage of 2M d ~ 1 vectors of size N d . In dimension 2 the classical method 
requires a storage of order 0(N 4 ) and our method requires a storage of order 0(MN 2 ). 
In dimension 3 the classical method requires a storage of order 0(N ) (thanks to the 
symmetries of the matrix of kernel modes, see |25| ) . and our method requires a storage of 
order 0(M 2 N 3 ). 

As a numerical example we report the results obtained in the case of space homogeneous 
two-dimensional Maxwellian molecules using as a comparison the exact analytic solution 
(see |41j). The results for the relative L\ norm of the error at time t = 0.01 are reported 
in Tabled 

Although further extensive testing is necessary, the results are very promising and seem 
to indicate a very low influence of the number of directions over the accuracy of the scheme. 
For M = 2 the angle error dominates, but as soon as M = 4 the error in N is dominating. 
Note that the number of angle directions will indirectly influence the aliasing effect trough 
the slight change in the relaxation times. This may explain the slight error variations that 
we observe taking M > 4. 
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Finally, in view of space non homogeneous computations, we will have the additional 
advantage of taking a larger number of gridpoints without increasing too much the com- 
putational cost, thus allowing the computations of flows at larger Mach number compared 
to conventional deterministic schemes. Further numerical results are under development 
and will be presented in the work |21j . 

4. Conclusions 

We have presented a deterministic way for computing the Boltzmann collision operator 
with fast algorithms, for a class of interactions which includes the case of hard spheres 
in dimension 3. The method is based on a Carleman-like representation of the operator 
that allows to express it as a combination of convolutions (this is trivially true for the 
loss part but it is not trivial for the gain part). A suitable periodized truncation of the 
operator is then used to derive new spectral methods computable with a high speed up 
in computation times. This brings the overall cost in dimension d to 0(M d ~ 1 N d log 2 N) 
where iV is the number of velocity parameters and M the number of angular directions 
in each angular coordinate. Consistency and accuracy of the proposed schemes are also 
presented, and it is shown to be spectrally accurate. Moreover the error on the momentum 
and energy is spectrally small and independent of the value of the speed-up parameter 
M. First numerical results seem to indicate the validity and the flexibility of the present 
approach that, to our opinion, will make deterministic schemes much more competitive 
with Monte Carlo methods in several situations. 

Appendix: Remarks on admissible collision kernels and an extension to 

the "non-decoupled" case 

Let us study the cases where the assumption (|3.2j) is satisfied. For hard spheres in 
dimension 3, or Maxwellian molecules in dimension 2, one has the equation ()3.2|) with 
a = b = 1. Formally for the Coulomb potential in dimension 3, we have 

B(0,\u\) = \u\~ 3 sin- 4 (6>/2), 

and thus, thanks to formula (|2.2jl 

B(x,y) = 2 d ~ 1 \x\~ 4 . 

This suggests, in dimension 3, to consider the following family of "variable hard sphere" 
collision kernels 

(0.4) B 7 (6, \u\) = sin 7 - 1 (0/2) \u\~< . 

Indeed simple computations give 

B 7 (x,y) =2 d - 1 |x| 7 - 1 

and thus they satisfy the decoupling assumption (|3,2j) . In the case where 7 G (—2, 1] the 
angular part of the collision kernel remains integrable. On the contrary, for 9 ~ 0, the 
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equivalent derived from the physical non explicit formula in ^1] for inverse-power laws 
kernels (for a potential with n such that 7 = (n — 5)/(n — 1)) is of the form 

B™ ct (e,\u\) ~ e ^ Q an 2 T L (9/2)\u\-> 

with 7 G [-3,1). It is therefore always non-integrable for ~ 0. 

The model (|0.4|) coincides with the hard spheres model for 7 = 1 and, formally, coincides 
with the kernel of the Coulomb potential for 7 = —3. Moreover for 7 6 (—2, 1] {i.e. hard 
potentials and the so-called moderately soft potentials) it remains integrable for 9 ~ 0. 
Thus it seems quite reasonable to consider it as a model for cutoff hard and moderately 
soft potentials, as well as hard spheres. 

In dimension 2 the same arguments and computations lead to the following cutoff hard 
and moderately soft potentials model 

B 7 (0,\u\) = sinT(0/2) \ u \^ 

valid for 7 £ (—3, 1], which coincides with the case of Maxwellian molecules for 7 = 0. 

For the spectral method the other situation where one obtains naturally a fast algorithm 
is the case where collisions concentrate on the grazing part: see |44| and |45j for a fast 
algorithm to compute the Fokker-Planck-Landau collision operator, which is the limit of 
the Boltzmann collision operator in the grazing collision limit. In this case indeed one 
of the two variables x or y of the representation (|2.4|) disappears in the limit process, 
which "decouples" the kernel modes. Thus it may be possible to construct fast algorithms 
for non-cutoff models by splitting the collision operator into a cutoff part treated by the 
method presented in this paper, and a non-cutoff part restricted to very small deviation 
angles, which would be close to the grazing collision limit and thus could be computed by 
the fast algorithm of gH 05] . 
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